The objective of this demo is to quantify the categories of herbaceous vegetation in San Diego city - USA. The dataset used is Copernicus Global Land Service (CGLS) which provides a series of bio-geophysical products on the status and evolution of land surface at global scale to a spatial resolution of 100m. More information click here

3.1 Requirements

library(rgee)
library(sf)
library(tidyverse)
library(mapview)
library(raster)
library(ggspatial)
ee_Initialize()

> ee_Initialize()
── rgee 1.1.2 ────────────────────────── earthengine-api 0.1.292 ── 
 ✓ user: not_defined
 ✓ Initializing Google Earth Engine:  DONE!
 ✓ Earth Engine account: users/ambarja 
───────────────────────────────────────────────────────────────────

3.2 Vector layer reading of your study area

sandiego <- st_read(
  "https://github.com/healthinnovation/sdb-gpkg/raw/main/SanDiego_districts.gpkg"
) %>%
  dplyr::select(PO_NAME) 
head(sandiego)

> head(sandiego)
Simple feature collection with 6 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -117.1162 ymin: 32.58161 xmax: -116.208 ymax: 32.90116
Geodetic CRS:  WGS 84
      PO_NAME                           geom
1      Alpine MULTIPOLYGON (((-116.6075 3...
2      Bonita MULTIPOLYGON (((-116.9655 3...
3   Boulevard MULTIPOLYGON (((-116.2085 3...
4       Campo MULTIPOLYGON (((-116.3573 3...
5 Chula Vista MULTIPOLYGON (((-117.0728 3...
6 Chula Vista MULTIPOLYGON (((-117.0081 3...

3.3 Cooking data with rgee

dataset <- ee$Image("COPERNICUS/Landcover/100m/Proba-V-C3/Global/2019")$
  select("discrete_classification")$
  eq(30) # 30: Herbaceous vegetation

herbaceous <- dataset$updateMask(dataset)

3.4 A function that extracts the number of pixels of herbaceous vegetation for each polygon

ee_get_data <- function(img, region, scale = 1000) {
  lista_histo <- list()
  for (i in 1:nrow(region)) {
    region_ee <- region[i, ] %>% sf_as_ee()
    ee_histo <- img$reduceRegion(
      reducer = ee$Reducer$frequencyHistogram(),
      geometry = region_ee,
      scale = scale
    )
    lista_histo[[i]] <- ee_histo$getInfo() %>%
      map_df(., .f = as.data.frame) %>%
      mutate(NAME = region[[i, 1]] %>% as.vector())
  }
  return(lista_histo)
}

3.5 Application of the function `ee_get_data`

# Time: ~ 2 min
rawdata <- ee_get_data(img = herbaceous, region = sandiego)
class(rawdata)

[1] "list"

3.6 Preprocesing rawdata

# Herbaceous area in Hectare Units 
newdata <- rawdata %>% 
  map_df(.f = as_tibble) %>% 
  mutate_if(
    is.numeric,
    .funs = function(x) {x * 1000 * 1000 / 10000}
    ) %>% 
  rename(area_Ha = X1)
head(newdata)

# A tibble: 6 × 2
  area_Ha NAME       
    <dbl> <chr>      
1   800   Alpine     
2  1073.  Bonita     
3   188.  Boulevard  
4  1300   Campo      
5    59.6 Chula Vista
6   179.  Chula Vista

# Left join with spatial data
sandiego_herb <- left_join(
 sandiego,
 newdata,
 by = c("PO_NAME"="NAME")
)
head(sandiego_herb)

Simple feature collection with 6 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -117.1162 ymin: 32.58161 xmax: -116.208 ymax: 32.90116
Geodetic CRS:  WGS 84
      PO_NAME    area_Ha                           geom
1      Alpine  800.00000 MULTIPOLYGON (((-116.6075 3...
2      Bonita 1073.33333 MULTIPOLYGON (((-116.9655 3...
3   Boulevard  187.84314 MULTIPOLYGON (((-116.2085 3...
4       Campo 1300.00000 MULTIPOLYGON (((-116.3573 3...
5 Chula Vista   59.60784 MULTIPOLYGON (((-117.0728 3...
6 Chula Vista  179.21569 MULTIPOLYGON (((-117.0728 3...

3.8 Exploration map

mapview(sandiego,zcol="area_Ha", layer.name = "Herbaceous")

3.9 Final Map

ggplot() + 
  geom_sf(data = sandiego_herb) +
  annotation_map_tile(zoom = 10) + 
  geom_sf(data = sandiego_herb,aes(fill = area_Ha),lwd=0.05) + 
  scale_fill_viridis() + 
  theme_minimal() + 
  labs(title = "Herbaceous vegetation - San Diego",
       caption = "Source: Copernicus Global Land Service (CGLS)")